d## How Economic Integration Affects Party Issue Emphases
## Graham, Kim, Tavits, Ward
## February 2015

## This file is to replicate 
## Graham, Kim, Tavits, and Ward ``How Economic Integration Affects Party Issue Emphases". 
## Comparative Political Studies, 48(10):1227-59, 2015.

# Set working directory
setwd("..")

## Load libraries
lapply(c("foreign", "MASS", "xtable", "plm", "texreg", "car", "arm", "lme4","plyr","apsrtable", "lmtest", "mvtnorm"), library, character.only=TRUE)

## Load data 
cmp <- read.dta("Data/Country_level.dta")
party <- read.dta("Data/Party_level.dta")
DD <- read.dta("Data/DD_data.dta")

############ Table 1 ############
## Model 1
table1_1 <- plm(unweighted_nonecon ~ kof_glob_econ + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced, model = "random", index = c("country_code", "election_count"), data = cmp)

## Model 2
table1_2 <- plm(unweighted_nonecon ~ EU_pc_econ_B + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced, model = "random", index = c("country_code", "election_count"), data = cmp)

## Model 3
table1_3 <- plm(econ_polar_lowe ~ kof_glob_econ + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced, model = "random", index = c("country_code", "election_count"), data = cmp)

## Model 4
table1_4 <- plm(econ_polar_lowe ~ EU_pc_econ_B + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced, model = "random", index = c("country_code", "election_count"), data = cmp)

############ Table 2 ############
## Model 1
table2_1 <- lmer(nonecon ~ kof_glob_econ + party_size + party_age + niche + incumbent_party + Religiosity2 + Ethnic_frac_Alesina + enep+ gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid  + advanced + (1|country_code) + (1|party), data=party)

## Model 2
table2_2 <- lmer(nonecon ~ EU_pc_econ_B + party_size + party_age + niche + incumbent_party  + Religiosity2 + Ethnic_frac_Alesina + enep+ gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced + (1|country_code) + (1|party), data=party)

############ Table 3: Difference-in-differences estimation ############

## Create a function that generates clustered standard errors in LM. 
clx <- function(fm, dfcw, cluster){
  library(sandwich)
  library(lmtest)
  M <- length(unique(cluster))
  N <- length(cluster)
  dfc <- (M/(M-1))*((N-1)/(N-fm$rank))
  u  <- apply(estfun(fm),2,function(x) tapply(x, cluster, sum))
  vcovCL <- dfc*sandwich(fm, meat=crossprod(u)/N)*dfcw
  coeftest(fm, vcovCL) }

## Model 1
table3_1 <- lm(nonecon ~ eu_member_ever + post + eu_member, data = DD)
clx(table3_1, 1, as.character( DD$party[as.numeric(rownames(table3_1$model))]))

## Model 2 
table3_2 <- lm(nonecon ~ eu_member_ever + post + eu_member + Ethnic_frac_Alesina + advanced + niche, data = DD)
clx(table3_2,1,as.character(DD$party[as.numeric(rownames(table3_2$model))]))

## Model 3
table3_3 <- lm( nonecon ~ eu_member_ever + post + eu_member + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + pop_pwt + gini_market_swiid + advanced + party_size + party_age + niche + incumbent_party, data= DD)
clx(table3_3,1,as.character(DD$party[as.numeric(rownames(table3_3$model))]))

## Model 4
table3_4 <- lm(nonecon ~ eu_member_ever * placebo_1999 + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + pop_pwt + gini_market_swiid + advanced + party_size + party_age + niche + incumbent_party, data=DD[DD$year  < 2004, ])
clx(table3_4,1,as.character(DD$party[as.numeric(rownames(table3_4$model))]))



############ Table SI.3.1 ############
## Create a function that generates clustered standard errors in PLM models. 
plmCluster <- function(mod_obj){
  G <- length(unique(unlist(attr(mod_obj[6]$model, "index")[1])))
  N <- length(unlist(attr(mod_obj[6]$model, "index")[1]))
  dfa <- (G / (G - 1)) * (N - 1)/mod_obj$df.residual
  return(dfa * vcovHC(mod_obj, type = "HC0", cluster = "group", adjust = TRUE))
}

## Create a function that returns the VCOV
clx2 <- function(fm, dfcw, cluster){
  library(sandwich)
  library(lmtest)
  M <- length(unique(cluster))
  N <- length(cluster)
  dfc <- (M/(M-1))*((N-1)/(N-fm$rank))
  u  <- apply(estfun(fm),2,function(x) tapply(x, cluster, sum))
  vcovCL <- dfc*sandwich(fm, meat=crossprod(u)/N)*dfcw
  return(vcovCL) }

## Model 1
tableSI3_1 <- plm(unweighted_nonecon ~ kof_glob_econ + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced, index = c("country_code", "election_count"), model = "random", data = cmp)
# Calculate clustered standard errors
plmCluster(tableSI3_1)

## Model 2
tableSI3_2 <- plm(unweighted_nonecon ~ EU_pc_econ_B + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced, model = "random", index = c("country_code", "election_count"), data = cmp)
# Calculate clustered standard errors
plmCluster(tableSI3_2)

## Model 3
tableSI3_3 <- plm(unweighted_nonecon ~ kof_glob_econ + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced, index = c("country_code", "election_count"), model = "within", data = cmp)
# Calculate clustered standard errors
plmCluster(tableSI3_3)

## Model 4
tableSI3_4 <- plm(unweighted_nonecon ~ EU_pc_econ_B + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced, model = "within", index = c("country_code", "election_count"), data = cmp)
# Calculate clustered standard errors
plmCluster(tableSI3_4)

## Model 5
tableSI3_5 <- lm(unweighted_nonecon ~  unweighted_nonecon_lag + kof_glob_econ + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced, data = cmp)
# Calculate clustered standard errors
round(clx(tableSI3_5, 1, as.character(cmp$countryname[as.numeric(rownames(tableSI3_5$model))])), digits = 3)

## Model 6
tableSI3_6 <- lm(unweighted_nonecon ~ unweighted_nonecon_lag + EU_pc_econ_B + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced, data = cmp)
# Calculate clustered standard errors
round(clx(tableSI3_6, 1, as.character(cmp$countryname[as.numeric(rownames(tableSI3_6$model))])), digits = 3)

############ Table SI.3.2 ############
## Model 1
tableSI3_7 <- plm(unweighted_nonecon ~ trade_wdi + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced, model = "random", index = c("country_code", "election_count"), data = cmp)
# Calculate clustered Standard Errors
tableSI3_7_SE <- sqrt(diag(tableSI3_7$vcov)) 

## Model 2
tableSI3_8 <- plm(unweighted_nonecon ~ trade_wdi + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced, model = "random", index = c("country_code", "election_count"), data = cmp)
# Calculate clustered standard errors
tableSI3_8_SE <- sqrt(diag(plmCluster(tableSI3_8)))

## Model 3
tableSI3_9 <- plm(unweighted_nonecon ~ trade_wdi + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced, model = "within", index = c("country_code", "election_count"), data = cmp)
# Calculate clustered standard errors
tableSI3_9_SE <- sqrt(diag(plmCluster(tableSI3_9)))

## Model 4
tableSI3_10 <- lm(unweighted_nonecon ~ unweighted_nonecon_lag + trade_wdi + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced, data = cmp)
# Calculate clustered standard errors
tableSI3_10_SE <- sqrt(diag(clx2(tableSI3_10, 1, as.character(cmp$countryname[as.numeric(rownames(tableSI3_10$model))]))))


############ Table SI.5.1 ############
party.data  <- ddply(party, .(election_count, party), summarize, nonecon = mean(nonecon, na.rm=T),
                     EU_pc_econ_B = mean(EU_pc_econ_B, na.rm=T), trade_wdi=mean(trade_wdi, na.rm=T),
                     kof_glob_econ = mean(kof_glob_econ, na.rm=T), party_size=mean(party_size, na.rm=T),
                     party_age=mean(party_age, na.rm=T), niche=mean(niche, na.rm=T), 
                     incumbent_party=mean(incumbent_party, na.rm=T),
                     Religiosity2=mean(Religiosity2, na.rm=T), 
                     Ethnic_frac_Alesina=mean(Ethnic_frac_Alesina, na.rm=T), enep=mean(enep, na.rm=T), 
                     gdp_growth_wdi=mean(gdp_growth_wdi, na.rm=T), pop_pwt=mean(pop_pwt, na.rm=T), 
                     gini_market_swiid=mean(gini_market_swiid, na.rm=T), advanced=mean(advanced, na.rm=T))
                     
## Model 1
tableSI5_1  <- plm(nonecon  ~ kof_glob_econ + party_size + party_age + niche + incumbent_party 
                + Religiosity2 + Ethnic_frac_Alesina + enep+ gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid 
                + advanced,  model = "random", index=c("party","election_count"), data=party.data)
coeftest(tableSI5_1, plmCluster(tableSI5_1)) # results reported in the table

## Model 2
tableSI5_2 <- plm(nonecon ~ EU_pc_econ_B + party_size + party_age + niche + incumbent_party 
                + Religiosity2 + Ethnic_frac_Alesina + enep+ gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid 
                + advanced,  model = "random", index = c("party", "election_count"), data=party.data)
coeftest(tableSI5_2, plmCluster(tableSI5_2)) # results reported in the table


## Model 3
tableSI5_3 <- plm(nonecon ~ kof_glob_econ + party_size + party_age + niche + incumbent_party 
                + Religiosity2 + Ethnic_frac_Alesina + enep+ gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid 
                + advanced,  model = "within", index = c("party", "election_count"), data=party.data)
coeftest(tableSI5_3, plmCluster(tableSI5_3)) # results reported in the table

## Model 4
tableSI5_4 <- plm(nonecon ~ EU_pc_econ_B + party_size + party_age + niche + incumbent_party 
                + Religiosity2 + Ethnic_frac_Alesina + enep+ gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid 
                + advanced,  model = "within", index = c("party", "election_count"), data=party.data)
coeftest(tableSI5_4, plmCluster(tableSI5_4)) # results reported in the table

## Model 5
tableSI5_5 <- lm(nonecon ~ lag_nonecon_PARTY + kof_glob_econ + party_size + party_age + niche + incumbent_party 
               + Religiosity2 + Ethnic_frac_Alesina + enep+ gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid 
               + advanced,  data=party)
# Calculate clustered SE
clx(tableSI5_5, 1, party$party[as.numeric(rownames(tableSI5_5$model))] )

## Model 6
tableSI5_6 <- lm(nonecon ~ lag_nonecon_PARTY + EU_pc_econ_B + party_size + party_age + niche + incumbent_party 
               + Religiosity2 + Ethnic_frac_Alesina + enep+ gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid 
               + advanced,  data=party)

# Calculate clustered SE
clx(tableSI5_6, 1, party$party[as.numeric(rownames(tableSI5_6$model))] )

############ Table SI.5.2 ############
# Model 1
tableSI5_7 <- lmer(nonecon ~ trade_wdi + party_size + party_age + niche + incumbent_party 
                 + Religiosity2 + Ethnic_frac_Alesina + enep+ gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid 
                 + advanced + (1|country_code) + (1|party), data=party)
# Model 2
tableSI5_8 <- plm(nonecon ~ trade_wdi + party_size + party_age + niche + incumbent_party 
                + Religiosity2 + Ethnic_frac_Alesina + enep+ gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid 
                + advanced,  model = "random", index = c("party", "election_count"), data=party.data)
coeftest(tableSI5_8, plmCluster(tableSI5_8)) # results reported in the table

# Model 3
tableSI5_9 <- plm(nonecon ~ trade_wdi + party_size + party_age + niche + incumbent_party 
                + Religiosity2 + Ethnic_frac_Alesina + enep+ gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid 
                + advanced,  model = "within", index = c("party", "election_count"), data=party.data)
coeftest(tableSI5_9, plmCluster(tableSI5_9)) # results reported in the table

# Model 4
tableSI5_10 <- lm(nonecon ~ lag_nonecon_PARTY + trade_wdi + party_size + party_age + niche + incumbent_party 
                + Religiosity2 + Ethnic_frac_Alesina + enep+ gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid 
                + advanced,  data=party)
# Clustered SE:
clx(tableSI5_10, 1, party$party[as.numeric(rownames(tableSI5_10$model))] )

############ Table SI.5.3 ############
## Model 1
tableSI5_11 <- lmer(nonecon ~ kof_glob_econ * niche + party_size + party_age + incumbent_party 
                 + Religiosity2 + Ethnic_frac_Alesina + enep+ gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid 
                 + advanced
                 + (1|country_code) + (1|party), data=party)

## Model 2
tableSI5_12 <- lmer(nonecon ~ kof_glob_econ * party_age + party_size + niche + incumbent_party 
                 + Religiosity2 + Ethnic_frac_Alesina + enep+ gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid 
                 + advanced
                 + (1|country_code) + (1|party), data=party)

## Model 3
tableSI5_13 <- lmer(nonecon ~ kof_glob_econ *party_size + party_age + niche + incumbent_party 
                 + Religiosity2 + Ethnic_frac_Alesina + enep+ gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid
                 + advanced
                 + (1|country_code) + (1|party), data=party)

## Model 4
tableSI5_144 <- lmer(nonecon ~ kof_glob_econ *incumbent_party + party_size + party_age + niche
                 + Religiosity2 + Ethnic_frac_Alesina + enep+ gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid 
                 + advanced
                 + (1|country_code) + (1|party), data=party)

####### Predicted outcomes for KOF model
#### Make a dataframe organized like the plm data!
countryelections <- attr(table1_1[6]$model, "index")
model <- as.data.frame(table1_1[6])
model$country_code <- countryelections$country_code
model$election_count <- countryelections$election_count

## use this to get the year variable we need into this data
model$id2 <- paste (as.character (model$country_code), as.character(model$election_count), sep="_") #Election yeard id variable
cmp$id2 <- paste(as.character (cmp$country_code), as.character(cmp$election_count), sep="_") #### new variable for making this data set
modelData <- merge(model, cmp[,c("id2", "countryname")], by = "id2", all.x = TRUE, sort = FALSE)
colnames(modelData) <- sapply(colnames(modelData), function(name) { gsub("model.","", name)})

#kof low
lowKOF <- data.frame(kof_glob_econ = mean(modelData$kof_glob_econ, na.rm = T) - sd(modelData$kof_glob_econ, na.rm = T),
                     Religiosity2 = mean(modelData$Religiosity2, na.rm = T),
                     Ethnic_frac_Alesina = mean(modelData$Ethnic_frac_Alesina, na.rm = T),
                     enep = mean(modelData$enep, na.rm = T), 
                     gdp_growth_wdi = mean(modelData$gdp_growth_wdi,na.rm = T),
                     gini_market_swiid = mean(modelData$gini_market_swiid, na.rm = T),
                     pop_pwt = mean(modelData$I.pop_pwt * 1000, na.rm = T),
                     advanced = 1)
#kof high
highKOF <- data.frame(kof_glob_econ = mean(modelData$kof_glob_econ, na.rm = T) + sd(modelData$kof_glob_econ, na.rm = T),
                      Religiosity2 = mean(modelData$Religiosity2, na.rm = T),
                      Ethnic_frac_Alesina = mean(modelData$Ethnic_frac_Alesina, na.rm = T),
                      enep = mean(modelData$enep, na.rm = T), 
                      gdp_growth_wdi = mean(modelData$gdp_growth_wdi,na.rm = T),
                      gini_market_swiid = mean(modelData$gini_market_swiid, na.rm = T),
                      pop_pwt = mean(modelData$I.pop_pwt * 1000, na.rm = T),
                      advanced = 1)

### Pg. 16: Predictions for the politicization model
predict(table1_1, newdata = lowKOF) #54.90871
predict(table1_1, newdata = highKOF) #60.88662

### Pg. 19: Predictions for the convergence/polarization model
predict(table1_3, newdata = lowKOF)  #2.199576
predict(table1_3, newdata = highKOF) #1.249414

####### Predicted outcomes for EU model

#### make a dataframe organized like the plm data!
countryelections <- attr(table1_2[6]$model, "index")
model <- as.data.frame(table1_2[6])
model$country_code <- countryelections$country_code
model$election_count <- countryelections$election_count

## use this to get the year variable we need into this data
model$id2 <- paste (as.character (model$country_code), as.character(model$election_count), sep="_") #Election yeard id variable
cmp$id2 <- paste(as.character (cmp$country_code), as.character(cmp$election_count), sep="_") #### new variable for making this data set
modelData <- merge(model, cmp[,c("id2", "countryname", "year")], by = "id2", all.x = TRUE, sort = FALSE)
colnames(modelData) <- sapply(colnames(modelData), function(name) { gsub("model.","", name)})

#EU non-member
lowEU <- data.frame(EU_pc_econ_B = min(modelData$EU_pc_econ_B),
                    Religiosity2 = mean(modelData$Religiosity2, na.rm = T),
                    Ethnic_frac_Alesina = mean(modelData$Ethnic_frac_Alesina, na.rm = T),
                    enep = mean(modelData$enep, na.rm = T), 
                    gdp_growth_wdi = mean(modelData$gdp_growth_wdi,na.rm = T),
                    gini_market_swiid = mean(modelData$gini_market_swiid, na.rm = T),
                    pop_pwt = mean(modelData$I.pop_pwt * 1000, na.rm = T),
                    advanced = 1)

#EU SEA
midEU <- data.frame(EU_pc_econ_B = 5.560110,
                    Religiosity2 = mean(modelData$Religiosity2, na.rm = T),
                    Ethnic_frac_Alesina = mean(modelData$Ethnic_frac_Alesina, na.rm = T),
                    enep = mean(modelData$enep, na.rm = T), 
                    gdp_growth_wdi = mean(modelData$gdp_growth_wdi,na.rm = T),
                    gini_market_swiid = mean(modelData$gini_market_swiid, na.rm = T),
                    pop_pwt = mean(modelData$I.pop_pwt * 1000, na.rm = T),
                    advanced = 1)

#EU most integrated
highEU <- data.frame(EU_pc_econ_B = max(modelData$EU_pc_econ_B),
                     Religiosity2 = mean(modelData$Religiosity2, na.rm = T),
                     Ethnic_frac_Alesina = mean(modelData$Ethnic_frac_Alesina, na.rm = T),
                     enep = mean(modelData$enep, na.rm = T), 
                     gdp_growth_wdi = mean(modelData$gdp_growth_wdi,na.rm = T),
                     gini_market_swiid = mean(modelData$gini_market_swiid, na.rm = T),
                     pop_pwt = mean(modelData$I.pop_pwt * 1000, na.rm = T),
                     advanced = 1)

# Pg. 17: Salience predictions
predict(table1_2, newdata = lowEU) # 56.16062
predict(table1_2, newdata = midEU) # 59.58151
predict(table1_2, newdata = highEU) # 61.61747


# Pg. 19: Convergence predictions
predict(table1_4, newdata = lowEU) # 1.914496
predict(table1_4, newdata = midEU) # 1.313995
predict(table1_4, newdata = highEU) # 0.9566042

############ Figure 1 ############
#use the same dataset generated in the above EU predictions section

#### calculate random effects 
alphaxbeta <- modelData[,2]-residuals(table1_2)
xbeta <- as.matrix(modelData[,3:10])%*%as.vector(coef(table1_2)[2:9])
random.intercepts <- alphaxbeta-xbeta
mean.deviations <- random.intercepts-coef(table1_2)[1]

#put them into this working data frame 
modelData$mean.dev <- mean.deviations
modelData$ran.int <- random.intercepts

#create a dataframe of just the 2004 ascension countries
expData <- modelData[modelData$countryname%in%c("Czech Republic", "Cyprus","Hungary","Malta","Latvia","Estonia","Slovenia","Lithuania","Slovakia","Poland"),]
expData[,c("countryname", "year")]
expData <- expData[as.character(c(305, 306, 310, 311, 318, 319, 327, 328, 340, 341, 357, 358, 364, 365)),]

#now make the predictive distributions 
set.seed(1801)
coef.sims <- mvrnorm(n = 1000, mu = c(mean(expData$ran.int), coef(table1_2)[2:9]), Sigma = vcov(table1_2))

estimate.y.pre <- vector()
estimate.y.post.no.eu <- vector()
estimate.y.post <- vector()

for(i in 1:1000){ 
  estimate.y.pre[i] <- coef.sims[i,]%*%c(1,apply(expData[c(2,4,6,8,10,12,14),3:10],2,mean))
  estimate.y.post.no.eu[i] <- coef.sims[i,]%*%c(1,min(expData$EU_pc_econ_B),apply(expData[c(1,3,5,7,9,11,13),4:10],2,mean))
  estimate.y.post[i] <-  coef.sims[i,]%*%c(1,apply(expData[c(1,3,5,7,9,11,13),3:10],2,mean))
}

opar <- par(mar=c(5,5,5,1), oma = c(1, 1, 1, 1))

plot(density(estimate.y.post),
     ylim=c(0,.24),
     xlim=c(43,68),
     yaxs="i",
     xaxs="i",
     xlab="Non-Economic Issues as a Percentage of Manifestos",
     main="",
     ylab="Kernal Density",
     xaxt="n",
     yaxt="n")
axis(1,at=c(45, 50, 55, 60, 65),labels=T,lwd=0,lwd.ticks=.5)
axis(2,at=c(.03,.06,.09,.12,.15,.18,.21),labels=T,lwd=0,lwd.ticks=.5)
lines(density(estimate.y.post.no.eu),lty=2)

lines(x=c(mean(estimate.y.post),mean(estimate.y.post)),y=c(0.09,.18))
lines(x=c(mean(estimate.y.post),mean(estimate.y.post)),y=c(0.01,.08))
text(x=c(mean(estimate.y.post),mean(estimate.y.post.no.eu)),y=c(.086,.066),labels=c("Mean","Mean"))
lines(x=c(mean(estimate.y.post.no.eu),mean(estimate.y.post.no.eu)),y=c(0.07,.22))
lines(x=c(mean(estimate.y.post.no.eu),mean(estimate.y.post.no.eu)),y=c(0.01,.06))

arrows(x1=mean(estimate.y.post),x0=mean(estimate.y.post),y1=0.002, y0=0.01,length=.15,lwd=1)
arrows(x1=mean(estimate.y.post),x0=mean(estimate.y.post),y1=0.002, y0=0.003,length=.15,lwd=1)
arrows(x1=mean(estimate.y.post),x0=mean(estimate.y.post),y1=0.002, y0=0.003,length=.15,angle=25,lwd=1)
arrows(x1=mean(estimate.y.post),x0=mean(estimate.y.post),y1=0.002, y0=0.003,length=.14,angle=20,lwd=1)
arrows(x1=mean(estimate.y.post),x0=mean(estimate.y.post),y1=0.002, y0=0.003,length=.13,angle=15,lwd=1)
arrows(x1=mean(estimate.y.post),x0=mean(estimate.y.post),y1=0.002, y0=0.003,length=.13,angle=10,lwd=1)
arrows(x1=mean(estimate.y.post),x0=mean(estimate.y.post),y1=0.002, y0=0.003,length=.13,angle=5,lwd=1)
arrows(x1=mean(estimate.y.post),x0=mean(estimate.y.post),y1=0.002, y0=0.003,length=.13,angle=4,lwd=1)
arrows(x1=mean(estimate.y.post),x0=mean(estimate.y.post),y1=0.002, y0=0.003,length=.13,angle=3,lwd=1)

arrows(x1=mean(estimate.y.post.no.eu),x0=mean(estimate.y.post.no.eu),y1=0.002, y0=0.01,length=.15,lwd=1)
arrows(x1=mean(estimate.y.post.no.eu),x0=mean(estimate.y.post.no.eu),y1=0.002, y0=0.003,length=.15,lwd=1)
arrows(x1=mean(estimate.y.post.no.eu),x0=mean(estimate.y.post.no.eu),y1=0.002, y0=0.003,length=.15,angle=25,lwd=1)
arrows(x1=mean(estimate.y.post.no.eu),x0=mean(estimate.y.post.no.eu),y1=0.002, y0=0.003,length=.14,angle=20,lwd=1)
arrows(x1=mean(estimate.y.post.no.eu),x0=mean(estimate.y.post.no.eu),y1=0.002, y0=0.003,length=.13,angle=15,lwd=1)
arrows(x1=mean(estimate.y.post.no.eu),x0=mean(estimate.y.post.no.eu),y1=0.002, y0=0.003,length=.13,angle=10,lwd=1)
arrows(x1=mean(estimate.y.post.no.eu),x0=mean(estimate.y.post.no.eu),y1=0.002, y0=0.003,length=.13,angle=5,lwd=1)
arrows(x1=mean(estimate.y.post.no.eu),x0=mean(estimate.y.post.no.eu),y1=0.002, y0=0.003,length=.13,angle=4,lwd=1)
arrows(x1=mean(estimate.y.post.no.eu),x0=mean(estimate.y.post.no.eu),y1=0.002, y0=0.003,length=.13,angle=2,lwd=1)

legend(x = 59, y=0.23,legend=c("No EU Membership","Nice-Era EU Membership"),lty=c(2,1),seg.len=1,bty="n",ncol=1,y.intersp=0.45)

############ Figure 2 ############
# Fit interaction models
int.mod1 <- lmer(nonecon ~ kof_glob_econ * niche + party_size + party_age + incumbent_party 
                 + Religiosity2 + Ethnic_frac_Alesina + enep+ gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid 
                 + advanced
                 + (1|country_code) + (1|party), data=party)

int.mod2 <- lmer(nonecon ~ kof_glob_econ * party_age + party_size + niche + incumbent_party 
                 + Religiosity2 + Ethnic_frac_Alesina + enep+ gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid 
                 + advanced
                 + (1|country_code) + (1|party), data=party)

int.mod3 <- lmer(nonecon ~ kof_glob_econ *party_size + party_age + niche + incumbent_party 
                 + Religiosity2 + Ethnic_frac_Alesina + enep+ gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid
                 + advanced
                 + (1|country_code) + (1|party), data=party)

int.mod4 <- lmer(nonecon ~ kof_glob_econ *incumbent_party + party_size + party_age + niche
                 + Religiosity2 + Ethnic_frac_Alesina + enep+ gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid 
                 + advanced
                 + (1|country_code) + (1|party), data=party)

# Niche=0
crit.t <- qt(0.975, df=2190-14)
lower0 <- fixef(int.mod1)[2] - crit.t*sqrt(summary(int.mod1)$vcov[2,2])
coef0 <- fixef(int.mod1)[2]
upper0 <- fixef(int.mod1)[2] + crit.t*sqrt(summary(int.mod1)$vcov[2,2])

# Niche=1
coef1 <- fixef(int.mod1)[2] + fixef(int.mod1)[14]
se1 <- sqrt(vcov(int.mod1)[2,2] + vcov(int.mod1)[14,14] + 2*vcov(int.mod1)[2,14]) 
lower1 <- coef1 - crit.t*se1
upper1 <- coef1 + crit.t*se1

# Upper-left plot:
par(mfrow=c(2,2), mar=c(4,4,1,0.5))
plot(1:5, xaxt="n", ylab="", type="n", ylim=c(-0.1,0.3),  xlab="")
axis(1, at=c(2,4), labels=c("Niche", "Mainstream"))
points(x=2, y=coef1, pch=20)
segments(2,lower1, 2, upper1)
points(x=4, y=coef0, pch=20)
abline(h=0, lty=2, col="grey")
segments(4,lower0, 4, upper0)
mtext("Marginal effect of globalization", side=2, cex=0.8, line=2)

# Incumbent=0
coef0 <- fixef(int.mod4)[2]
lower0 <-coef0 - crit.t*sqrt(vcov(int.mod4)[2,2])
upper0 <-coef0 + crit.t*sqrt(vcov(int.mod4)[2,2])
# Incumbent=1
coef1 <- fixef(int.mod4)[2] + fixef(int.mod4)[14]
se1 <- sqrt(vcov(int.mod4)[2,2] + vcov(int.mod4)[14,14] + 2*vcov(int.mod4)[2,14])
lower1 <- coef1 - crit.t*se1
upper1 <- coef1 + crit.t*se1

# Upper-right plot:
plot(1:5, xaxt="n", ylab="", type="n", ylim=c(-0.1,0.3),  xlab="")
axis(1, at=c(2,4), labels=c("Incumbent", "Opposition"))
points(x=2, y=coef1, pch=20)
segments(2,lower1, 2, upper1)
points(x=4, y=coef0, pch=20)
segments(4,lower0, 4, upper0)
abline(h=0, lty=2, col="grey")
mtext("Marginal effect of globalization", side=2, cex=0.8, line=2)

# Bottom-left plot:
age.quantile <- quantile(party.data$party_age, na.rm=TRUE)
ruler <- seq(0, age.quantile[4], by=0.1)
effects <- fixef(int.mod2)[2] + fixef(int.mod2)[14]*ruler
se <- sqrt(vcov(int.mod2)[2,2] + ruler^2*vcov(int.mod2)[14,14] + 2*ruler*vcov(int.mod2)[2,14])

plot(x=ruler, y=effects, type="l", xlab="", ylab="",  ylim=c(-0.1,0.4))
lines(x=ruler, y=effects-crit.t*se, lty=2, col="grey")
lines(x=ruler, y=effects+crit.t*se, lty=2, col="grey")
mtext("Party Age", cex=0.8, line=2, side=1)
mtext("Marginal effect of globalization", side=2, cex=0.8, line=2)
legend("topleft", lty=2, col="grey", bty="n", legend="95% CI")

# Bottom-right plot:
size.quantile <- quantile(party.data$party_size, na.rm=TRUE)
ruler <- seq(0, 70, by=0.1)
effects <- fixef(int.mod3)[2] + fixef(int.mod3)[14]*ruler
se <- sqrt(vcov(int.mod3)[2,2] + ruler^2*vcov(int.mod3)[14,14] + 2*ruler*vcov(int.mod3)[2,14])
plot(x=ruler, y=effects, type="l", xlab="", ylab="",  ylim=c(-0.1,0.4))
lines(x=ruler, y=effects-crit.t*se, lty=2, col="grey")
lines(x=ruler, y=effects+crit.t*se, lty=2, col="grey")
mtext("Party Size", cex=0.8, line=2, side=1)
mtext("Marginal effect of globalization", side=2, cex=0.8, line=2)
legend("topleft", lty=2, col="grey", bty="n", legend="95% CI")


######## SI pg. 4: 1) Excluding Governmental and Administrative Efficiency item in DV
SI_pg_4_1 <- plm(unweighted_nonecon_no_gv_ffcncy ~ kof_glob_econ + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced, model = "random", index = c("country_code", "election_count"), data = cmp)

######## SI pg. 4: 2) Excluding Education expansion/limitation item in DV
SI_pg_4_2 <- plm(unweighted_nonecon_no_education ~ kof_glob_econ + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced, model = "random", index = c("country_code", "election_count"), data = cmp)

######## SI pg. 4: 3) Excluding culture item in DV
SI_pg_4_3 <- plm(unweighted_nonecon_no_culture ~ kof_glob_econ + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced, model = "random", index = c("country_code", "election_count"), data = cmp)

######## Footnote 10: Using "unweighted_nonecon_cee" as DV
FN_10_1 <- plm(unweighted_nonecon_cee ~ kof_glob_econ + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced, model = "random", index = c("country_code", "election_count"), data = cmp)
FN_10_2 <- plm(unweighted_nonecon_cee ~ EU_pc_econ_B + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced, model = "random", index = c("country_code", "election_count"), data = cmp)

######## Footnote 12-1: Controlling for manifesto length
FN_12_1 <- plm(unweighted_nonecon ~ kof_glob_econ + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced + total, model = "random", index = c("country_code", "election_count"), data = cmp)

######## Footnote 12-2: Using "relative non_econ" as DV
FN_12_2 <- plm(relative_nonecon ~ kof_glob_econ + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced + total, model = "random", index = c("country_code", "election_count"), data = cmp)

######## Footnote 13: Using weighted non_econ as DV
FN_13 <- plm(weighted_nonecon ~ kof_glob_econ + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced + total, model = "random", index = c("country_code", "election_count"), data = cmp)

######## Footnote 14: Using alternativemeasure of globalization - trade
FN_14_1 <- plm(unweighted_nonecon ~ trade_wdi + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced + total, model = "random", index = c("country_code", "election_count"), data = cmp)
FN_14_2 <- plm(unweighted_nonecon ~ openc_pwt + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced + total, model = "random", index = c("country_code", "election_count"), data = cmp)
FN_14_3 <- plm(unweighted_nonecon ~ openk_pwt + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced + total, model = "random", index = c("country_code", "election_count"), data = cmp)

######## Footnote 15: Using alternativemeasure of globalization - KOF overall measure
FN_15 <- plm(unweighted_nonecon ~ kof_glob_overall + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced + total, model = "random", index = c("country_code", "election_count"), data = cmp)

######## Footnote 16: Using alternativemeasure of European integration
FN_15_1 <- plm(unweighted_nonecon ~ EU_pc_all_B + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced + total, model = "random", index = c("country_code", "election_count"), data = cmp)
FN_15_2 <- plm(unweighted_nonecon ~ eu_pc_all_hm + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced + total, model = "random", index = c("country_code", "election_count"), data = cmp)
FN_15_3 <- plm(unweighted_nonecon ~ eu_jolly + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced + total, model = "random", index = c("country_code", "election_count"), data = cmp)

######## Footnote 19: Using DV excluding EU-related items
FN_19 <- plm(unweighted_nonecon_no_EU ~ kof_glob_econ + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_market_swiid + advanced + total, model = "random", index = c("country_code", "election_count"), data = cmp)

######## Footnote 23: DD Analaysis, no Bulgaria and Romania
FN_23 <- lm( nonecon ~ eu_member_ever + post + eu_member + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + pop_pwt + gini_market_swiid + advanced + party_size + party_age + niche + incumbent_party, data= DD[ !DD$countryname %in% c("Bulgaria", "Romania"),])
clx(FN_23,1,as.character(DD$party[as.numeric(rownames(FN_23$model))]))

######## Footnote SI-1: gini net measure in the place of gini gross
FN_SI_1 <- plm(unweighted_nonecon ~ kof_glob_econ + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_net_swiid + advanced, model = "random", index = c("country_code", "election_count"), data = cmp)

######## Footnote SI-3: Alternative religiosity measures
FN_SI_3_1 <- plm(unweighted_nonecon ~ kof_glob_econ + religiosity_wvs + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_net_swiid + advanced, model = "random", index = c("country_code", "election_count"), data = cmp)
FN_SI_3_2 <- plm(unweighted_nonecon ~ kof_glob_econ + ChurchAttend2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_net_swiid + advanced, model = "random", index = c("country_code", "election_count"), data = cmp)

######## Footnote SI-4: Controlling for inflation
FN_SI_4 <- plm(unweighted_nonecon ~ kof_glob_econ + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_net_swiid + advanced + inflation_conpri_wdi, model = "random", index = c("country_code", "election_count"), data = cmp)

######## Footnote SI-5: Controlling for democracy using Polity IV data
FN_SI_5 <- plm(unweighted_nonecon ~ kof_glob_econ + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + I(pop_pwt/1000) + gini_net_swiid + polity2, model = "random", index = c("country_code", "election_count"), data = cmp)

######## Footnote SI-7: DD Analysis, only countries with elections in both periods
FN_SI_7 <- lm( nonecon ~ eu_member_ever + post + eu_member + Religiosity2 + Ethnic_frac_Alesina + enep + gdp_growth_wdi + pop_pwt + gini_market_swiid + advanced + party_size + party_age + niche + incumbent_party, data= DD[DD$no_2004_elections == 0, ])
clx(FN_SI_7, 1, as.character(DD$party[ as.numeric(rownames(FN_SI_7$model)) ]))

